%{
AUTHOR: Felipe Arteaga

DATE:  2020--
-------------------------------------------------------------------------
PROJECT:
-------------------------------------------------------------------------
DESCRIPTION:
[Analsis de solo 2018 y 2019, 2017 se hace aparte y 2020 no hay una
población que recibe SMS muy grande]

Analiza la población SMS, que son los que postularon mientras el pop-up no
estaba activo y no volvieron a postular con el pop-up activo al menos hasta
el momento que se definió enviar los SMSs. Entre el 2017 y 2020 son los
siguientes:




2017:
SMS para no evaluados. 4 versiones randomizadas, diferentes cutoff
(analisis en analisisRD_SMS_2017.m)

2018:
SMSA para no evaluados. SMSB al final para todavía riesgosos (que no
recibieron SMSA)


2019:
SMS para no evaluados. SMS1 para no RM, SMS2 para RM, enviados en distintos
días. No se envió recordatorio

2020:
Solo se envió recordatorio (dos veces, una Martín desde Mienduc y otra
nosotros, que fue una mezcla entre el experimento whatsapp y SMSs). La
población SMS son los que no fueron evaluados por pop-up (grados no de
entrada de la RM), pero por ahora no los estoy considerando.

=========================================================================
%}

clc;clear;close all;fclose('all');
feature('DefaultCharacterSet', 'UTF-8');

pcName=char(java.lang.System.getProperty('user.name'));
if(strcmp(pcName,'felipe'))
    % PC Felipe
    myDir='/Users/felipe/Dropbox/';
    projectDir=[myDir,'git/warnings/'];
    projectDirData=[myDir,'projects/warnings/'];
end


addpath(genpath([myDir,'/myMatlabFunctions/']));

%%
dirPlots=[projectDir,'/paper/figuresCL/RDs/'];

% warning('Alternative dir plots')
% dirPlots=[projectDir,'/paper/figuresCL/RDs/Aux/'];

dirPlotsR=[projectDir,'/paper/figuresCL/RDs_multBandwidth/'];
dirTable=[projectDir,'/paper/tablesCL/'];

dirData=[projectDirData,'/data/chile/'];

heter=''; %'' , 'vulnerable', 'maturity '.  '': none;
closePlotsWhilePlotting=true;

plotRDs=false;
savePlotRDs=false;
plotWithZoom=false;
bandwidthTable='local'; % 'full','local','localAuto';


% For making subfigures with plots:
maxHorzPlots=3;
maxVertPlots=4;

tableBandwidthComparison=false; % Make a table comparing estimates with different bandwidths.
posBandwidth=1; % This is for the bandwidht robust excercise. Indicates the subpopulation that we are using the bandwidth table
plotAndSaveRobust=false; % Plots of robust bandwidht RD

tableDetails=false;
posDetails=1;% This is for the details RD table. Indicates the subpopulation that we are using for the details table

if(not(isempty(heter)))
    heter=['_',heter];
    plotRDs=false;
    tableBandwidthComparison=false;
    plotAndSaveRobust=false;
    tableDetails=false;
end
switch bandwidthTable
    case 'full'
        bwnote='Full bandwidth was used (-0.28,+0.68) with 2nd order (local) polynomial fit. ';
    case 'local'
        bwnote='Bandwidth was set to +-0.1 with default local polynomial fit (1st order). ';
    case 'localAuto'
        bwnote='Bandwidth was set to automatic (bwselect with default options) with default polynomial fit (1st order). ';
end

%% Load data:
anho=0; %0: pooled
if(anho==0)
    anhoStr='';
else
    anhoStr=sprintf('%i',anho);
end


notIn2020={}; % {'assignedToPref','acceptOffer','declineOffer','preferenciaAsign_1','placedInAnyPreferenceIniEnd','placedInAddedIniEnd','placedInAddedNoWaitlistIniEnd','placedInOldIniEnd','enrolledInAssigned','enrolledInAddedIniEnd','enrolledInAddedNoWaitlistIniEnd','enrolledInAddedNoRiskIniEnd'};

if(anho==0)
    
    d1=load([dirData,'2018/inputRD'],'dataRD');
    d2=load([dirData,'2019/inputRD'],'dataRD');
    d3=load([dirData,'2020/inputRD'],'dataRD');
    % Variables that are not available for 2020 yet:
    fillWithNan=notIn2020;
    %fillWithNan=
    for v=1:length(fillWithNan)
        d3.dataRD.(fillWithNan{v})=nan(height(d3.dataRD),1);
    end
    
    
    % Universe that is comparable between 2018-2020
    entryGrades=[-1 0 1 7 9];
    sirven1=ismember(d1.dataRD.grade,entryGrades)&not(d1.dataRD.cod_reg==13);
    sirven2=ismember(d2.dataRD.grade,entryGrades)&not(d2.dataRD.cod_reg==13);
    sirven3=ismember(d3.dataRD.grade,entryGrades)&not(d3.dataRD.cod_reg==13);
    
    d1.dataRD.comparable=sirven1;
    d2.dataRD.comparable=sirven2;
    d3.dataRD.comparable=sirven3;
    
    d1.dataRD.anho=2018*ones(height(d1.dataRD),1);
    d2.dataRD.anho=2019*ones(height(d2.dataRD),1);
    d3.dataRD.anho=2020*ones(height(d3.dataRD),1);
    
    
    % Genero riskSMS
    
    d1.dataRD.riskSMS=d1.dataRD.riskSMSA;
    d1.dataRD.pobSMS=d1.dataRD.pobSMSA;
    d1.dataRD.treatedSMS=d1.dataRD.treatedSMSA;
    
    d2.dataRD.riskSMS=d2.dataRD.riskSMS1;
    d2.dataRD.riskSMS(d2.dataRD.pobSMS2)=d2.dataRD.riskSMS2(d2.dataRD.pobSMS2);
    d2.dataRD.pobSMS=d2.dataRD.pobSMS1|d2.dataRD.pobSMS2;
    d2.dataRD.treatedSMS=d2.dataRD.treatedSMS1;
    d2.dataRD.treatedSMS(d2.dataRD.pobSMS2)=d2.dataRD.treatedSMS2(d2.dataRD.pobSMS2);
    
    
    
     vars={'changeAny','addAny','schoolsAdded','DRiskExpost','placedInAnyPreference'};
     for v=1:length(vars)
         d1.dataRD.([vars{v},'SMS'])=d1.dataRD.([vars{v},'PreSAPreSB']);
         d2.dataRD.([vars{v},'SMS'])=d2.dataRD.([vars{v},'PreS1End']);
         d2.dataRD.([vars{v},'SMS'])(d2.dataRD.pobSMS2)=d2.dataRD.([vars{v},'PreS2End'])(d2.dataRD.pobSMS2);
     end
     
     
     
     
    % Keep vars in common
    
    %incommon=intersect(intersect(d1.dataRD.Properties.VariableNames,d2.dataRD.Properties.VariableNames),d3.dataRD.Properties.VariableNames);
    incommon=intersect(d1.dataRD.Properties.VariableNames,d2.dataRD.Properties.VariableNames);
    incommon=incommon(not(strcmp(incommon,'mrun')));
    %dataRD=[d1.dataRD(:,incommon);d2.dataRD(:,incommon);d3.dataRD(:,incommon)];
    dataRD=[d1.dataRD(:,incommon);d2.dataRD(:,incommon)];
    
    
else
    load([dirData,anhoStr,'/inputRD'],'dataRD')
end

dataRD.rural=dataRD.newMarket<0;

% As asignment at .3 is control, just move (if any) in .3 to epsilon to the
% left:




dataRD.riskSMS(dataRD.riskSMS==.3)=.29999999;

% Avoid mass points at extremes:
dataRD.restrAll=dataRD.pobSMS==1&dataRD.riskSMS>.01&dataRD.riskSMS<.99;

% No need to be more specific, and put "Predicted placement risk of 1st
% attempt"
dataRD.Properties.VariableDescriptions{'riskSMS'}='Predicted placement risk';
dataRD.Properties.VariableDescriptions{'treatedSMS'}='Received SMS';



   
   
   


%% RD

% Defines subpopulation in which I want to calculate RDs for outcomes
% defined soon.

if(anho==0) % This calculates estimates pooling all years
    
    %     subpobs=cell(1,3);
    %     s=1;
    %     subpobs(s,:)={'2018',   dataRD.restrAll&dataRD.anho==2018,'2018'};s=s+1;
    %     subpobs(s,:)={'2019',      dataRD.restrAll&dataRD.anho==2019,'2019'};s=s+1;
    %     subpobs(s,:)={'2020',    dataRD.restrAll&dataRD.anho==2020,'2020'};s=s+1;
    %     subpobs(s,:)={'AddMandatory',            dataRD.restrAll&dataRD.appLength1Erased==1,'obligadosAAgregar'};s=s+1;
    %     subpobs(s,:)={'NotAddMandatory',            dataRD.restrAll&dataRD.appLength1Erased==0,'noOligadosAAgregar'};s=s+1;
    %     subpobs(s,:)={'All',            dataRD.restrAll,'pooled'};s=s+1;
    
    %     subpobs=cell(1,3);
    %     s=1;
    %     subpobs=cell(1,3);
    %     subpobs(s,:)={'2018 AddMand',   dataRD.restrAll&dataRD.anho==2018&dataRD.appLength1Erased==1,'2018addMand'};s=s+1;
    %     subpobs(s,:)={'2019 AddMand',      dataRD.restrAll&dataRD.anho==2019&dataRD.appLength1Erased==1,'2019addMand'};s=s+1;
    %     subpobs(s,:)={'2020 AddMand',    dataRD.restrAll&dataRD.anho==2020&dataRD.appLength1Erased==1,'2020addMand'};s=s+1;
    %     subpobs(s,:)={'2018 NotAddMand',   dataRD.restrAll&dataRD.anho==2018&dataRD.appLength1Erased==0,'2018notAddMand'};s=s+1;
    %     subpobs(s,:)={'2019 NotAddMand',      dataRD.restrAll&dataRD.anho==2019&dataRD.appLength1Erased==0,'2019notAddMand'};s=s+1;
    %     subpobs(s,:)={'2020 NotAddMand',    dataRD.restrAll&dataRD.anho==2020&dataRD.appLength1Erased==0,'2020notAddMand'};s=s+1;
    %     subpobs(s,:)={'All AddMand',            dataRD.restrAll&dataRD.appLength1Erased==1,'pooledAddMand'};s=s+1;
    %     subpobs(s,:)={'All NotAddMand',            dataRD.restrAll&dataRD.appLength1Erased==0,'poolednotAddMand'};s=s+1;
    
    % 4th column is if we calculate IV
    
    switch heter
        case ''
            subpobs=cell(1,4);
            s=1;
            subpobs(s,:)={'All',    dataRD.restrAll,'pooled','addAnySMS'};s=s+1;
            subpobs(s,:)={'2018',   dataRD.restrAll&dataRD.anho==2018,'2018',''};s=s+1;
            subpobs(s,:)={'2019',      dataRD.restrAll&dataRD.anho==2019,'2019',''};s=s+1;
            %subpobs(s,:)={'2020',    dataRD.restrAll&dataRD.anho==2020,'2020',''};s=s+1;
            posAll=1;
        
        otherwise
            error('aca')
    end
    
   
    
else
    subpobs=cell(1,3);
    s=1;
    subpobs(s,:)={'All',            dataRD.restrAll,'all',''};s=s+1;
    
    %     subpobs(2,:)={'Vulnerable',   dataRD.restrAll&dataRD.esSep==1,'sep'};
    %     subpobs(3,:)={'Non-vulnerable',      dataRD.restrAll&dataRD.esSep==0,'sepNo'};
    %     subpobs(4,:)={'PK',    dataRD.restrAll&dataRD.grade==-1,'gPk'};
    %     subpobs(5,:)={'9th',dataRD.restrAll&dataRD.grade==9,'g9'};
    %     subpobs(6,:)={'Voluntary',   dataRD.restrAll&dataRD.voluntary==1,'vol'};
    %     subpobs(7,:)={'Non-voluntary',      dataRD.restrAll&dataRD.voluntary==0,'volNo'};
end

assert(posBandwidth==posAll)
assert(posDetails==posAll)

% Definition of outcomes for RDs. 3rd row defines the beginning of a panel
% in the table
% 4th row


% New: enrolled conditional on placed (short name because of stata)
dataRD.enrolledInAssignedCOA=double(dataRD.enrolledInAssigned);
dataRD.enrolledInAssignedCOA(dataRD.assignedToPref==0)=nan;
dataRD.placedInAnyPreferenceNWIE=dataRD.placedInAnyPreferenceNoWaitlistIniEnd;
dataRD.placedInAnyPreferenceNoRiskIE=dataRD.placedInAnyPreferenceNoRiskIniEnd;

%notIn2020=[notIn2020,'enrolledInAssignedCOA'];

% 5th column defines if IV is calculated or not.
outcomes={'esSep','Economically Vulnerable','sep','A. Balance',0;...
    'rural','Rural','rur','',0;...
    'changeAnySMS','Any modification','cha','B. Choice Behavior',0;...
    'addAnySMS','Add any','aa','',0;...
    'schoolsAddedSMS','Schools Added','sa','',1;...
    % 'riskLastDayIni','True Risk initial attempt','ri','';...
    % 'riskLastDayEnd','True Risk final attempt','rf','';...
    % 'DRiskLastDayIniEnd','$\Delta$ True Risk','dr','';...
    %'riskExpostIni','True Risk initial attempt','riEP','';...
    %'riskExpostEnd','True Risk final attempt','rfEF','';...
    'DRiskExpostSMS','$\Delta$ Risk','drEP','',1;...
    %'numberOfNewUncongestedIniEnd','Uncongested schools Added','saUncong','';...
    %'addAsFirstIniEnd','Add as first','af','',1;...
    %'addInBetweenIniEnd','Add to middle','ab','',1;...
    %'addAsLastIniEnd','Add as last','al','',1;...
    %'deleteAnyIniEnd','Drop any','da','',1;...
    %'changeOrigOrderIniEnd','Re-order','ro','',1;...
    % 'riskLastDayEnd','Risk final portfolio','rf','';...
    'placedInAnyPreferenceSMS','Placed to preference','ap','C. Choice outcome',1;...
    %'placedInAddedIniEnd','Placed to added preference','apad','';...
    %'placedInOldIniEnd','Placed to old preference','apold','';...
    'enrolledInAssigned','Enrolled in placed','ea','',1;...
    %'enrolledInAssignedCOA','Enrolled in placed$|$placed','eaadcp','',1;...
    %'distanceEnrolled','Distance to enrolled (km)','d','',1;...
    %'valueAddedEnrolled','Value added enrolled','va','',1;...
    %'declineOffer','Decline offer','do','',1;...
    %'enrolledInAddedIniEnd','Enrolled in added','eaad','';...
    %'addAnyUncongestedIniEnd','Add any uncongested','aaUncong','D. Congestion-related outcomes',1;...
    %    'deltaProbPlacedNoWaitlistsIniEnd','$\Delta$ prob. placed to uncongested','dpw','',1;...
     %'deltaProbPlacedNoRiskIniEnd','$\Delta$ prob. placed to zero risk','dpnr','',1;...
    %'placedInAnyPreferenceNWIE','Placed to uncong.','apUncong','',1;...
    %'placedInAddedNoWaitlistIniEnd','Placed to uncong. added pref.','apadUncong','',1;...
    %'enrolledInAddedNoWaitlistIniEnd','Enrolled in uncong. added','eaadUncong','',1;...
%         'addAnyNoRiskIniEnd','Add any zero risk','aaNR','D2. Risk-related outcomes',1;...
%     'placedInAnyPreferenceNoRiskIE','Placed to zero risk','apNR','',1;...
%     'placedInAddedNoRiskIniEnd','Placed to zero risk added pref.','apadNR','',1;...
%     'enrolledInAddedNoRiskIniEnd','Enrolled in zero risk added','eaadNR','',1;...
    %'acceptOffer','Accept offer','ao','';...
    %'preferenciaAsign_1','Order assigned','oa','';...
    };


% outcomes={   % 'riskLastDayIni','True Risk initial attempt','ri','';...
%     %'riskLastDayEnd','True Risk final attempt','rf','';...
%     %'DRiskLastDayIniEnd','$\Delta$ True Risk','dr','';...
%     'riskExpostIni','True Risk initial attempt','riEP','';...
%     'riskExpostEnd','True Risk final attempt','rfEF','';...
%     'DRiskExpostIniEnd','$\Delta$ True Risk ExP','drEP','';...
%     'placedInAnyPreferenceIniEnd','Placed to preference','ap','C. Choice outcome';...
%     'placedInAddedIniEnd','Placed to added preference','apad','';...
%     'placedInOldIniEnd','Placed to old preference','apold','';...
%     };


if(anho==2020)
    outcomes=outcomes(not(ismember(outcomes(:,1),notIn2020)),:);
end



dataRD.Properties.VariableDescriptions{'riskSMS'}='Predicted placement risk';
for o=1:size(outcomes,1)
    
    dataRD.Properties.VariableDescriptions{outcomes{o,1}}=outcomes{o,2};
end


assert(allunique(outcomes(:,1)))
assert(allunique(outcomes(:,2)))
assert(allunique(outcomes(:,3)))


data=dataRD(:,[{'riskSMS'},outcomes(:,1)']);
inAnyPob=false(height(data),1);


% Generate stata commands
commands=cell(size(subpobs,1)*size(outcomes,1),2);
counter=1;

for p=1:size(subpobs,1)
    
    data.(sprintf('subpop%i',p))=subpobs{p,2};
    inAnyPob=inAnyPob|subpobs{p,2};
    
    for o=1:size(outcomes,1)
        
        % Avoid computing outcomes that are not available:
        noCalcular=(strcmp(subpobs{p,3},'2020')|strcmp(subpobs{p,3},'2020comp'))&ismember(outcomes{o,1},notIn2020);
        
        if(not(noCalcular))
            
            
            % Full bandwidth
            if(tableBandwidthComparison||strcmp(bandwidthTable,'full')||plotRDs)
                commands(counter,:)={sprintf('%s_%i_full',outcomes{o,3},p),sprintf('rdrobust %s riskSMS if subpop%i==1, c(.3) p(2) h(.28 .68)',outcomes{o,1},p)};counter=counter+1;
            end
            % Full loccal fixed bandwidth
            if(tableBandwidthComparison||strcmp(bandwidthTable,'local'))
                commands(counter,:)={sprintf('%s_%i_local',outcomes{o,3},p),sprintf('rdrobust %s riskSMS if subpop%i==1, c(.3)  h(.1 .1)',outcomes{o,1},p)};counter=counter+1;
            end
            % Full local auto bandwidht
            if(tableBandwidthComparison||strcmp(bandwidthTable,'localAuto'))
                commands(counter,:)={sprintf('%s_%i_localAuto',outcomes{o,3},p),sprintf('rdrobust %s riskSMS if subpop%i==1, c(.3) ',outcomes{o,1},p)};counter=counter+1;
            end
            
            % If with IV:
            if(outcomes{o,5}==1&&~isempty(subpobs{p,4}))
                
                endogenousVar=subpobs{p,4};
                
                % Full bandwidth
                if(strcmp(bandwidthTable,'full'))
                    commands(counter,:)={sprintf('%s_%i_full_iv',outcomes{o,3},p),sprintf('rdrobust %s riskSMS if subpop%i==1, fuzzy(%s) c(.3) p(2) h(.28 .68)',outcomes{o,1},p,endogenousVar)};counter=counter+1;
                end
                % Full loccal fixed bandwidth
                if(strcmp(bandwidthTable,'local'))
                    commands(counter,:)={sprintf('%s_%i_local_iv',outcomes{o,3},p),sprintf('rdrobust %s riskSMS if subpop%i==1, fuzzy(%s) c(.3)  h(.1 .1)',outcomes{o,1},p,endogenousVar)};counter=counter+1;
                end
                % Full local auto bandwidht
                if(strcmp(bandwidthTable,'localAuto'))
                    commands(counter,:)={sprintf('%s_%i_localAuto_iv',outcomes{o,3},p),sprintf('rdrobust %s riskSMS if subpop%i==1, fuzzy(%s) c(.3) ',outcomes{o,1},p,endogenousVar)};counter=counter+1;
                end
            end
        end
    end    
end
commands=commands(not(cellfun(@isempty,commands(:,1),'UniformOutput',true)),:);


if(false)
    % Run specific rd
dataRD.auxRest=dataRD.restrAll&dataRD.mandatoryToAdd==0;
a=stataCommand('rdrobust lengthEnd riskSMS if auxRest==1, c(.3) p(1) h(.1 .1)',dataRD)
stataExpress('rdrobust  DRiskExpostIniEnd riskSMS if auxRest==1, fuzzy(addAnyIniEnd) c(.3) p(1) h(.1 .1)',data)
end

% Run Stata:
data=data(inAnyPob,:);
res=stataCommand(commands,data);


%% Make plots and table


cantPobs=(size(subpobs,1));
cantRegs=(size(outcomes,1));

matrix_nl=nan(cantRegs,cantPobs);
matrix_nr=nan(cantRegs,cantPobs);
matrix_b=nan(cantRegs,cantPobs);
matrix_se=nan(cantRegs,cantPobs);

matrix_b_iv=nan(cantRegs,cantPobs);
matrix_se_iv=nan(cantRegs,cantPobs);
matrix_nl_iv=nan(cantRegs,cantPobs);
matrix_nr_iv=nan(cantRegs,cantPobs);

matrix_b_l=nan(cantRegs,cantPobs);
matrix_se_l=nan(cantRegs,cantPobs);

matrix_b_r=nan(cantRegs,cantPobs);
matrix_se_r=nan(cantRegs,cantPobs);

matrix_biasCorr_ci=nan(cantRegs,cantPobs,2);



regNames=fieldnames(res);

for p=1:cantPobs
    for r=1:cantRegs
        
        regName_i=sprintf('%s_%i_%s',outcomes{r,3},p,bandwidthTable);
        if(ismember(regName_i,regNames))
            res_i=res.(regName_i);
            
            % Beta
            matrix_b(r,p)=res_i.tau_cl;
            % Sd Beta
            matrix_se(r,p)=res_i.se_tau_cl;
            
            % Beta left
            matrix_b_l(r,p)=res_i.beta_p_l(1);
            % Sd Beta left
            matrix_se_l(r,p)=sqrt(res_i.V_cl_l(1));
            
            % Beta right
            matrix_b_r(r,p)=res_i.beta_p_r(1);
            % Sd Beta right
            matrix_se_r(r,p)=sqrt(res_i.V_cl_r(1));
            
            % Bias corrected Confidence interval (alpha=5%)
            matrix_biasCorr_ci(r,p,1)=res_i.tau_bc-norminv(.975)*res_i.se_tau_rb;
            matrix_biasCorr_ci(r,p,2)=res_i.tau_bc+norminv(.975)*res_i.se_tau_rb;
            
            % N
            matrix_nl(r,p)=res_i.N_h_l;
            matrix_nr(r,p)=res_i.N_h_r;
            
            % If with IV:
            if(outcomes{r,5}==1&&~isempty(subpobs{p,4}))
                regName_i_iv=sprintf('%s_%i_%s_iv',outcomes{r,3},p,bandwidthTable);
                
                res_i_iv=res.(regName_i_iv);
                % Beta
                matrix_b_iv(r,p)=res_i_iv.tau_cl;
                % Sd Beta
                matrix_se_iv(r,p)=res_i_iv.se_tau_cl;
                
                % N
                matrix_nl_iv(r,p)=res_i_iv.N_h_l;
                matrix_nr_iv(r,p)=res_i_iv.N_h_r;
                
            end
            
            if(plotRDs)
                
                % Load full bandwidth result for ploting
                warning('Add small bandwidth estimate for full bandwidth plot')
                
                res_i=res.(sprintf('%s_%i_%s',outcomes{r,3},p,'full'));
                res_i_table=res.(sprintf('%s_%i_%s',outcomes{r,3},p,bandwidthTable));
                
                plotsub=subpobs{p,2};
                
                figure
                plotRD(res_i,dataRD,plotsub,'otherPointEstimate',res_i_table);
                if(savePlotRDs)
                    easyExport([dirPlots,sprintf(sprintf('%s_S_%s_%s.png',anhoStr,outcomes{r,3},subpobs{p,3}))]);
                end
                
                
                
                if(plotWithZoom)
                    % Plot with zoom
                    figure
                    plotRD(res_i,dataRD,plotsub,'newxlim',[.1,.5],'nb',100,'otherPointEstimate',res_i_table);
                    easyExport([dirPlots,sprintf(sprintf('%s_S_%s_%s_withZoom.png',anhoStr,outcomes{r,3},subpobs{p,3}))]);
                end
                
                
                if(closePlotsWhilePlotting)
                    close all
                end
            else
                fprintf('Ojo: reg %s no existe para subpob %s\n',regName_i,subpobs{p,1});
            end
        end
    end
    

    
end





%% RD main table
matTable=nan(size(matrix_b).*[1 2]+[2 0]);
matTable_sd=nan(size(matrix_b).*[1 2]+[2 0]);

matTable(:,1:2:end)=[matrix_b;...
    matrix_nl(end,:)...
    ;matrix_nr(end,:)];

matTable(:,2:2:end)=[matrix_b_iv;...
    matrix_nl_iv(find(cell2mat(outcomes(:,5)),1,'first'),:)...
    ;matrix_nr_iv(find(cell2mat(outcomes(:,5)),1,'first'),:)];


matTable_sd(:,1:2:end)=[matrix_se;nan(2,size(matrix_se,2))];
matTable_sd(:,2:2:end)=[matrix_se_iv;nan(2,size(matrix_se_iv,2))];

withInfo=not(all(isnan(matTable),1));


header=repmat({''},2,size(matTable,2));
header(1,1:2:end)=subpobs(:,1)';
header(1,2:2:end)=subpobs(:,1)';
header(2,2:2:end)={'IV'};

matTable=matTable(:,withInfo);
matTable_sd=matTable_sd(:,withInfo);
header=header(:,withInfo);


opt=struct;
opt.header=header;
opt.stderrs=mat2cellstr(matTable_sd,'conParentesis',true,'precisionDecimal','%.3f');
opt.primeracolumna=[outcomes(:,2);{'NL';'NR'}];
%opt.stars=getStars(matTable,matTable_sd);
opt.addColumnNumber=true;
if(anho==0)
    opt.columnafantasma=4;

else
    %opt.columnafantasma=[4];
    opt.title=sprintf('%s Sample - SMS effect',anhoStr);
end
opt.filaFantasma=size(matTable,1)-2;
paneles=find(not(ismissing(outcomes(:,end-1))));
opt.panel=[num2cell(paneles-1),outcomes(paneles,end-1)];
%opt.alignmentFirstCol={'L{3cm}'};
opt.adjust=true;

opt.file=sprintf('%s%s_tableSMS%s',dirTable,anhoStr,heter);


opt.title='RD estimates of SMSs effects';
opt.label='tab:RD-SMS';

opt.note=['Local linear RD estimates of warning''s SMS effects. Computed using triangular kernel with bandwidth 0.1. Heteroskedasticity-robust nearest neighbor variance estimator with minimum of 3 neighbors reported in parentheses; computed as in \citet{calonico2014robust}. ',...
    'We report estimates in the pooled sample and for each year. IV (column 2) shows the instrumental variable specifications, where the endogenous regressor is the add any school indicator. ',...
    'Panel A: predetermined covariates. Panel B: measures of choice behavior from initial to final application. $\Delta$ risk is change in application risk from first to final attempt. ``Add to X'''' are additions of schools in given place on list, relative to initial application submission. ',...
    'Panel C: outcomes of choice process. ``Enrolled in placed'''' is equal to one for students who receive a placement and enroll in the placed school. ``Enrolled in placed | placed'''' is the enrollment rate in the placed school, conditional on receiving a placement. ',...
    'Panel D: congestion attributes of behavior and placement outcomes. ``Uncongested'''' schools are those with excess capacity. '];


tabla=cell2latex(mat2cellstr(matTable,'precisionDecimal','%.3f','revisarFilas',true),'opts',opt);
compileLatex(tabla)



%% RD details table
% Only pooled resutls
% Columns:
% 1) Estimate
% 4-5) Bias-corrected confidence interval (95\%)
% 2) Estimate left
% 3) Estimate right

if(tableDetails)
ci=cellfun(@(x,y)['[',x,' ',y,']'],...
    mat2cellstr(matrix_biasCorr_ci(:,posDetails,1),'precisionDecimal','%.3f'),...
    mat2cellstr(matrix_biasCorr_ci(:,posDetails,2),'precisionDecimal','%.3f'),...
    'uniformOutput',false);
    

cellTable=[mat2cellstr(matrix_b(:,posDetails),'precisionDecimal','%.3f'),...
    ci,...
    mat2cellstr(matrix_b_l(:,posDetails),'precisionDecimal','%.3f'),...
    mat2cellstr(matrix_b_r(:,posDetails),'precisionDecimal','%.3f'),...
    mat2cellstr(100*(matrix_b_r(:,posDetails)-matrix_b_l(:,posDetails))./matrix_b_l(:,posDetails),'precisionDecimal','%.1f','sufijo','%')];
    
matTable_sd=[matrix_se(:,posDetails),nan(size(matrix_se,1),1),matrix_se_l(:,posDetails),matrix_se_r(:,posDetails),nan(size(matrix_se,1),1)];


opt=struct;
opt.header={'Estimate','Bias-corrected C.I.','Estimates','Estimates','%';...
        '','','Left','Right','change'};
opt.stderrs=mat2cellstr(matTable_sd,'conParentesis',true,'precisionDecimal','%.3f');
opt.primeracolumna=outcomes(:,2);
%opt.stars=getStars(matTable,matTable_sd);
opt.addColumnNumber=true;
if(anho==0)
    opt.columnafantasma=2;
    opt.title='Details of RD Pooled Pop-Up effects';
else
    %opt.columnafantasma=[4];
    opt.title=sprintf('%s Sample - Pop-Up effect',anhoStr);
end

paneles=find(not(ismissing(outcomes(:,4))));
opt.panel=[num2cell(paneles-1),outcomes(paneles,4)];
%opt.alignmentFirstCol={'L{3cm}'};
opt.adjust=true;
opt.label='tab:RD-SMS-details';
opt.file=sprintf('%s%s_tableSMS_details%s',dirTable,anhoStr,heter);
opt.note=['Local linear RD estimates of warning''s SMS effects. Computed using triangular kernel with bandwidth 0.1. Heteroskedasticity-robust nearest neighbor variance estimator with minimum of 3 neighbors reported in parentheses; computed as in \citet{calonico2014robust}. ',...
    'We report estimate in the pooled sample and splitting by year. ',...
    'Panel A: predetermined covariates. Panel B: measures of choice behavior from initial to final application. $\Delta$ risk is change in application risk from first to final attempt. ``Add to X'''' are additions of schools in given place on list, relative to initial application submission. ',...
    'Panel C: measures outcomes of choice process. ``Enrolled in placed'''' is equal to one for students who receive a placement and enroll in the placed school. ``Enrolled in placed | placed'''' is the enrollment rate in the placed school, conditional on receiving a placement. ',...
    'Panel D: congestion attributes of behavior and placement outcomes. ``Uncongested'' schools are those with excess capacity. '];
tablaD=cell2latex(cellTable,'opts',opt);
compileLatex(tablaD)
end

%% Make table of bandwidth comparison

if(tableBandwidthComparison)
    cantRobuts=3;
    labelRobust={'Full','+-0.1','rdbwselect'};
    cantRegs=(size(outcomes,1));
    
    matrix_nl=nan(cantRegs,cantRobuts);
    matrix_nr=nan(cantRegs,cantRobuts);
    matrix_b=nan(cantRegs,cantRobuts);
    matrix_se=nan(cantRegs,cantRobuts);
    
    matrix_bl=nan(cantRegs,cantRobuts);
    matrix_br=nan(cantRegs,cantRobuts);
    
    regNames=fieldnames(res);
    
    
    cellsubfig=cell(maxHorzPlots,ceil(cantRegs/maxHorzPlots)); % Defined transposed to fill it with scalar indexing.
    counter=1;
    for p=posBandwidth
        for r=1:cantRegs
            
            regName_i=sprintf('%s_%i_full',outcomes{r,3},p);
            regName_i_local=sprintf('%s_%i_local',outcomes{r,3},p);
            regName_i_localAuto=sprintf('%s_%i_localAuto',outcomes{r,3},p);
            
            if(ismember(regName_i_local,regNames))
                res_i=cell(cantRobuts,1);
                
                res_i{1}=res.(regName_i);
                res_i{2}=res.(regName_i_local);
                res_i{3}=res.(regName_i_localAuto); % AUTO tiene que ser la ultima!
                
                
                for i=1:cantRobuts
                    % Beta
                    matrix_b(r,i)=res_i{i}.tau_cl;
                    % Sd Beta
                    matrix_se(r,i)=res_i{i}.se_tau_cl;
                    
                    % N
                    matrix_nl(r,i)=res_i{i}.N_h_l;
                    matrix_nr(r,i)=res_i{i}.N_h_r;
                    
                    % N
                    matrix_bl(r,i)=res_i{i}.h_l;
                    matrix_br(r,i)=res_i{i}.h_r;
                    
                end
                
                if(plotAndSaveRobust)
                    
                    plotsub=subpobs{p,2};
                    
                    
                    
                    % Plot alternatives bandwidths:
                    
                    colors=linspecer(3);
                    c1=colors(1,:);
                    c2=colors(2,:);
                    c3=colors(3,:);
                    
                    figure
                    % Plot full bandwidth:
                    
                    a1=plotRD(res_i{1},dataRD,plotsub,'color',c1); hold on;
                    if(a1.posBeta>.5);shift=-.1;else;shift=.1;end
                    a2=plotRD(res_i{2},dataRD,plotsub,'withCutoffVerticalLine',0,'withBinsreg',0,'posbeta',a1.posBeta+shift,'color',c2);  hold on;
                    a3=plotRD(res_i{3},dataRD,plotsub,'withCutoffVerticalLine',0,'withBinsreg',0,'posbeta',a1.posBeta+2*shift,'color',c3);
                    
                    legend([a1.functionLine,a2.functionLine,a3.functionLine],labelRobust,'interpreter','latex')
                    
                    hold off
                    
                    % Guardo pal latexSubfigure
                    cellsubfig{counter}=easyExport([dirPlotsR,sprintf(sprintf('%s_S_%s_%s_bandComparison%s.png',anhoStr,outcomes{r,3},subpobs{p,3}))],'caption',outcomes{r,2});
                    counter=counter+1;
                    
                    
                    close
                    
                    %                 % Plot with zoom
                    %                 figure
                    %                 a=plotRD(res_i,dataRD,plotsub,'newxlim',[.1,.5]);
                    %                 easyExport([dirPlots,sprintf(sprintf('%s_S_%s_%s_withZoom.png',anhoStr,outcomes{r,3},subpobs{p,3}))]);
                    %
                    %                 close
                    
                end
            end
        end
        
        
    end
    
   %%
    if(plotAndSaveRobust)
    for p=1:ceil(size(cellsubfig,2)/maxVertPlots)
        
        newEnd=min((p*maxVertPlots),size(cellsubfig,2));
        note='SMS warning RD effect fits and point estimates by bandwidth for outcomes listed in panel titles. ``Full'''': global quadratic. ``+/- 0.1'''': local linear within 0.1 bandwidth. ``rdbwselect'''': optimal bandwidth selection using \citet{calonico2014robust}. See  section \ref{sec:RD} for details.';
        if p==1
            label='fig:appBW';
        else
            label=sprintf('fig:appBW_%i',p);
        end
        subfiguresLatex(cellsubfig(:,((p-1)*maxVertPlots+1):newEnd)','file',[dirPlotsR,'subfigRobustRDs_',num2str(p)],'erp','figuresCL/RDs_multBandwidth/','erpf', 'figuresCL/RDs_multBandwidth/','caption','Multiple bandwidths RD plots of platform-based pop-up warning effects',...
            'scale',.52,'docwidth',1.2,'note',note,'label',label);
    end
    end
    %%
    
    cellTable=[[mat2cellstr(matrix_b,'precisionDecimal','%.3f','revisarCols',true),mat2cellstr([matrix_bl(:,end),matrix_br(:,end),matrix_nl(:,end),matrix_nr(:,end)],'precisionDecimal','%.2f','revisarCols',true)];...
        mat2cellstr([[matrix_nl(1,1:2);matrix_nr(1,1:2)],nan(2,5)])];
    
    nans=nan(size(cellTable));
    nans(1:size(matrix_se,1),1:size(matrix_se,2))=matrix_se;
    matTable_sd=nans;
    
    
    opt=struct;
    opt.header=[[{'Bandwidth'},labelRobust,{'rdbwselect','rdbwselect','rdbwselect','rdbwselect'}];...
        {'','Estimate','Estimate','Estimate','BW left','BW right','N left','N right'}];
    opt.stderrs=mat2cellstr(matTable_sd,'conParentesis',true,'precisionDecimal','%.3f');
    opt.primeracolumna=[outcomes(:,2);{'N left';'N right'}];
    %opt.stars=getStars(matTable,matTable_sd);
    opt.addColumnNumber=true;
    opt.columnafantasma=[1 2];
    opt.title='Platform SMS RD estimates with alternate bandwidths';
    opt.filaFantasma=size(outcomes,1);
    paneles=find(not(ismissing(outcomes(:,4))));
    opt.panel=[num2cell(paneles-1),outcomes(paneles,4)];
    %opt.alignmentFirstCol={'L{3cm}'};
    opt.adjust=true;
    opt.label='tab:SMS-BW';
    opt.file=sprintf('%s%s_tableSMS-Bandwidth%s',dirTable,anhoStr,heter);
    opt.note=['Local linear and full sample quadratic polynomial RD estimates of warning''s SMS effects.. Computed using triangular kernel with different bandwidths. ``Full'''' bandwidth uses 2nd order polynomial fit, while "+-0.1" and rdbwselect uses 1st order (local) polynomial. Heteroskedasticity-robust nearest neighbor variance estimator with minimum of 3 neighbors reported in parentheses; computed as in \citet{calonico2014robust}. ',...
    'Panel A: predetermined covariates. Panel B: measures of choice behavior from initial to final application. $\Delta$ risk is change in application risk from first to final attempt. ``Add to X'''' are additions of schools in given place on list, relative to initial application submission. ',...
    'Panel C: measures outcomes of choice process. ``Enrolled in placed'''' is equal to one for students who receive a placement and enroll in the placed school. ``Enrolled in placed | placed'''' is the enrollment rate in the placed school, conditional on receiving a placement. ',...
    'Panel D: congestion attributes of behavior and placement outcomes. ``Uncongested'' schools are those with excess capacity. '];
    tablaComp=cell2latex(cellTable,'opts',opt);
    compileLatex(tablaComp)
    
end

% %% Sync dropbox with Overleaf folder
% if(strcmp(pcName,'felipe'))
% overleafDir='/Users/felipe/Dropbox/Apps/Overleaf/Warnings'' draft/';
% gitDir='/Users/felipe/Dropbox/git/warningsChris/Presentations/Slides_AEA_Jan2021/FiguresSlides/';
% dirPlotsOverleaf=[gitDir,'/figuresRDsChi/'];
% dirPlotsROverleaf=[gitDir,'/figuresRDsChi_robust/'];
% dirTableOverleaf=[overleafDir,'/tablesChi/'];
%
% system(sprintf('rsync -avzh "%s" "%s" --exclude=".*"',dirTable,dirTableOverleaf));
% system(sprintf('rsync -avzh "%s" "%s" --exclude=".*"',dirPlots,dirPlotsOverleaf));
% system(sprintf('rsync -avzh "%s" "%s" --exclude=".*"',dirPlotsR,dirPlotsROverleaf));
%
% end